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Abstract  A  coupled  ocean-atmosphere  mesoscale  ensem¬ 
ble  prediction  system  has  been  developed  by  the  Naval 
Research  Laboratory.  This  paper  describes  the  components 
and  implementation  of  the  system  and  presents  baseline 
results  from  coupled  ensemble  simulations  for  two  tropical 
cyclones.  The  system  is  designed  to  take  into  account  major 
sources  of  uncertainty  in:  (1)  non-deterministic  dynamics, 
(2)  model  error,  and  (3)  initial  states.  The  purpose  of  the 
system  is  to  provide  mesoscale  ensemble  forecasts  for  use 
in  probabilistic  products,  such  as  reliability  and  frequency 
of  occurrence,  and  in  risk  management  applications.  The 
system  components  include  COAMPS®  (Coupled  Ocean/ 
Atmosphere  Mesoscale  Prediction  System)  and  NCOM 
(Navy  Coastal  Ocean  Model)  for  atmosphere  and  ocean 
forecasting  and  NAVDAS  (NRL  Atmospheric  Variational 
Data  Assimilation  System)  and  NCODA  (Navy  Coupled 
Ocean  Data  Assimilation)  for  atmosphere  and  ocean  data 
assimilation.  NAVDAS  and  NCODA  are  3D-variational 
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(3D VAR)  analysis  schemes.  The  ensembles  are  generated 
using  separate  applications  of  the  Ensemble  Transform  (ET) 
technique  in  both  the  atmosphere  (for  moving  or  non¬ 
moving  nests)  and  the  ocean.  The  atmospheric  ET  is 
computed  using  wind,  temperature,  and  moisture  variables, 
while  the  oceanographic  ET  is  derived  from  ocean  current, 
temperature,  and  salinity  variables.  Estimates  of  analysis 
error  covariance,  which  is  used  as  a  constraint  in  the  ET,  are 
provided  by  the  ocean  and  atmosphere  3DVAR  assimilation 
systems.  The  newly  developed  system  has  been  successfully 
tested  for  a  variety  of  configurations,  including  differing 
model  resolution,  number  of  members,  forecast  length,  and 
moving  and  fixed  nest  options.  Results  from  relatively  coarse 
resolution  (~27-km)  ensemble  simulations  of  Hurricanes 
Hanna  and  Ike  demonstrate  that  the  ensemble  can  provide 
valuable  uncertainty  information  about  the  storm  track  and 
intensity,  though  the  ensemble  mean  provides  only  a  small 
amount  of  improved  predictive  skill  compared  to  the 
deterministic  control  member. 

Keywords  Ensembles  •  Coupled  ocean-atmosphere 
modeling  •  Tropical  cyclones  •  Variational  analyses  • 
Probabilistic  prediction 

1  Introduction 

All  models  are  imperfect  and  all  forecasts  are  inaccurate. 
Hence,  at  best  we  can  hope  to  forecast  a  distribution  of 
possible  future  states  given  our  model  and  a  history  of 
observations.  A  compelling  method  for  attempting  to 
sample  this  distribution  is  ensemble  forecasting  (Toth  and 
Kalnay  1993).  The  idea  is  to  make  an  ensemble  of  forecasts 
whose  differing  initial  conditions  describe  our  best  estimate 
of  the  uncertainty  of  the  initial  state  estimate  and  whose 
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differing  representations  of  sub-grid  scale  processes  repre¬ 
sent  our  uncertainty  in  the  model.  Rigorously  accurate 
ensemble  methods  for  describing  the  distribution  of  future 
states  given  past  information  include  particle  filters  and 
Monte  Carlo-Markov  chains  (Van  Leeuwen  2009).  Unfor¬ 
tunately,  the  computational  resources  required  for  the 
application  of  such  methods  to  models  with  tens  of  millions 
of  variables  such  as  the  coupled  ocean-atmosphere  system 
considered  are  not  readily  available.  For  this  reason, 
computationally  efficient  approximations  to  these  schemes 
must  be  used.  Here  we  employ  a  framework  in  which  the 
mean  of  the  ensemble  of  initial  conditions  for  the  coupled 
ocean-atmosphere  state  is  obtained  from  three-dimensional 
variational  (3DVAR)  schemes  in  the  atmosphere  and  ocean 
and  the  ensemble  is  obtained  by  inputting  estimates  of 
analysis  error  variance  from  the  3D  VAR  data  assimilation 
schemes  into  the  Ensemble  Transform  (ET)  ensemble 
generation  technique  of  Bishop  and  Toth  (1999),  McLay 
et  al.  (2008),  and  Bishop  et  al.  (2009). 

The  advent  of  regional  two-way  coupled  models  is 
relatively  recent  and  perhaps  because  of  this  very  little 
work  has  been  done  in  the  area  of  regional  coupled  model 
ensemble  forecasting.  It  is,  nevertheless,  a  very  promising 
area  of  research  because  uncertainty  in  forecasts  of  the 
oceanic  boundary  layer  is  directly  linked  to  uncertainty  in 
forecasts  of  the  atmospheric  boundary  layer  (ABL)  and  vice 
versa.  It  is  extremely  difficult  to  model  this  interaction 
between  uncertainty  in  the  ocean  and  atmosphere  without  a 
coupled  model.  This  paper  represents  one  of  the  first 
attempts  to  describe  this  process. 

A  description  of  the  components  of  the  coupled 
ensemble  system  is  given  in  “Section  2”.  “Section  3” 
describes  the  ET  system,  and  “Section  4”  details  the 
implementation  of  the  mesoscale  ensemble  system.  Some 
basic  results  for  two  tropical  cyclone  (TC)  case  studies  are 
provided  in  “Section  5”  to  provide  a  baseline  for  the 
ensemble  system.  We  choose  TC  case  studies  because  of  the 
strong  air-ocean  coupling  in  TCs  and  the  need  for 
improvement  in  both  track  and  intensity  forecasts.  While 
TC  track  forecasts  have  improved  over  the  past  30  years  in 
the  Atlantic  basin,  with  average  track  forecast  errors 
reduced  noticeably  (Rappaport  et  al.  2009),  it  remains  a 
critical  issue  to  increase  the  accuracy  of  TC  track  forecasts 
in  the  days  prior  to  potential  landfall  to  provide  lead  time 
for  better  protection  of  coastal  communities.  In  the  2008 
hurricane  season,  Hanna  and  Ike  were  two  consecutive 
storms  that  made  landfall  on  the  east  and  southeast  coasts 
of  the  USA,  respectively.  The  coupled  ensemble  system  is 
used  to  simulate  these  two  storms  and  make  an  initial 
assessment  of  system  performance  for  both  track  and 
intensity  forecasts  over  a  2-  to  3 -day  period  before  landfall. 
The  paper  concludes  with  a  summary  and  conclusions 
(“Section  6”). 


2  Description  of  the  ocean-atmosphere  coupled 
ensemble  prediction  system 

The  coupled  system  developed  at  the  Naval  Research 
Laboratory  (NRL)  is  comprised  of  four  primary  components, 
as  described  in  Chen  et  al.  (2010):  the  nonhydrostatic 
atmospheric  model  COAMPS  (Coupled  Ocean/Atmosphere 
Mesoscale  Prediction  System;  Hodur  1997),  the  hydrostatic 
ocean  model  NCOM  (Navy  Coastal  Ocean  Model;  Martin 
2000),  the  atmospheric  data  assimilation  system  NAVDAS 
(NRL  Variational  Data  Assimilation  System;  Daley  and  Barker 
2001a,  b),  and  the  ocean  data  assimilation  system  NCODA 
(Navy  Coupled  Ocean  Data  Assimilation  System;  Cummings 
2005).  A  detailed  description  of  the  components  is  given  in 
these  references,  so  only  a  brief  overview  with  emphasis  on 
the  new  modifications  to  the  system  relevant  to  the  ensemble 
implementation  applied  to  TC  applications  is  provided  here. 

2.1  COAMPS 

The  physical  parameterizations  in  the  atmospheric  component 
include  a  cumulus  parameterization  (used  for  grid  spacing 
greater  than  10  km)  based  on  a  modified  Kain-Fritsch  (K-F) 
(Kain  and  Fritsch  1993;  Kain  2004)  scheme,  a  cloud 
microphysics  scheme  based  on  Rutledge  and  Hobbs  (1983), 
and  a  radiation  parameterization  based  on  Harshvardhan  et  al. 
(1987).  The  sea  surface  temperature  (SST)  feedback  is 
through  a  Monin-Obukhov-based  surface  similarity  scheme. 
The  surface  momentum,  heat,  and  moisture  fluxes  are 
parameterized  using  standard  bulk  aerodynamic  formulae 
(Louis  1979): 


Ui  =  CqU2, 

(i) 

u*e*  =  CnUAe, 

(2) 

U*q*  =  C^U  Aq , 

(3) 

where  U*  is  the  friction  velocity,  0 *  is  the  convective 
scaling  temperature,  q*  is  the  mixing  ratio  scaling,  U  is  the 
mean  surface  wind  speed,  AO  and  Aq  are  the  mean 
potential  temperature  and  vapor  mixing  ratio  differences 
computed  using  the  surface  and  the  first  model  grid  point 
above  the  surface  (~10  m),  and  CD,  CH,  and  CE  are  the  bulk 
transfer  coefficients  for  momentum,  heat,  and  moisture, 
respectively.  The  transfer  coefficients  are  stability  and 
roughness  length  dependent,  as  derived  from  observations, 
with  differing  momentum  roughness  length  as  compared  to 
heat  and  moisture.  The  stability  function  is  a  modified 
version  of  Louis  (1979)  based  on  version  3.0  of  the 
Tropical  Ocean  and  Global  Atmosphere  Coupled  Ocean- 
Atmosphere  Response  Experiment  (TOGA  CO  ARE;  Fairall 
et  al.  2003). 
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Recent  enhancements  to  several  of  the  CO  AMPS 
physical  parameterizations  that  reflect  an  increased 
understanding  from  recent  TC  observations  and  research 
have  been  made.  The  changes  include  a  reduction  of  CD 
under  high  winds,  the  inclusion  of  heating  due  to 
dissipation  of  turbulent  kinetic  energy  (TKE),  and  the 
decrease  of  ice  nuclei  concentration  to  reduce  upper-level 
cloud  ice.  Laboratory  tank  experiments  suggest  that  CD 
asymptotically  approaches  2.5  xlCT3  when  the  surface 
wind  speed  exceeds  35  m  s-1  (Donelan  et  al.  2004).  This 
phenomenon  has  also  been  observed  in  wind  profile 
measurements  in  TCs  (Powell  et  al.  2003)  and  direct 
turbulence  measurement  at  near-hurricane-force  wind 
speeds  during  the  latest  Coupled  Boundary  Layers  Air- 
Sea  Transfer  (CBLAST)  hurricane  field  experiment, 
suggesting  an  even  lower  asymptotic  value  of  23  m  s-1 
(Black  et  al.  2007).  However,  the  CBLAST  measurements 
displayed  a  large  fluctuation  of  CD  values  for  a  given  1 0-m 
wind  speed  and  had  no  CD  values  for  winds  exceeding 
30  m  s_1.  Therefore,  the  neutral  10-m  CD  over  the  ocean  in 
COAMPS  has  been  modified  to  be  in  agreement  with  that 
estimated  by  Donelan  et  al.  (2004).  Another  change  is 
related  to  heating  caused  by  energy  loss  due  to  TKE 
dissipation.  This  process,  parameterized  and  tested  in 
COAMPS,  shows  a  positive  impact  on  TC  intensity  and 
structure  forecasts  (Jin  et  al.  2007).  More  recently,  the  ice 
nucleation  formulation  has  been  modified  based  on 
previous  studies.  This  change  results  in  better  TC  structure 
as  well  as  intensity  forecasts. 

2.2  NCOM 


The  hydrostatic,  mesoscale  version  of  the  Navy  Coastal 
Ocean  Model  (NCOM)  with  hybrid  sigma-z  levels  is  the 
ocean  component  of  the  coupled  system.  Lorcing  fields 
from  the  atmospheric  component  are  used  as  ocean  surface 
boundary  conditions  for  momentum,  potential  temperature, 
and  salinity  as  shown  by  Eqs.  4-7  (Martin  2000): 


Po’ 


(4) 


Ku 


dv 

dz 


Ty_ 

Po’ 


(5) 


<90  2b  +  Qq  +  Qs 

Kn^~  = - 1  (6) 

dz  p0cp 

Kk-  =  S(Ew-Pt),  (7) 


where  KM  and  Ku  are  vertical  eddy  coefficients  for  the 
momentum  and  scalar  fields,  respectively;  u  and  v  are  the 
horizontal  components  of  the  current  fields,  0  and  S  are  the 
potential  temperature  and  salinity,  respectively;  rx  and  ry  are 
the  surface  wind  stress  in  the  x  and  y  directions, 
respectively;  Qh ,  Qe ,  and  Qs  are  the  net  long  wave,  latent, 
and  sensible  surface  heat  fluxes,  respectively;  Ey  and  PT  are 
the  surface  evaporation  and  precipitation  rates,  respectively; 
and  po  and  cp  are  the  density  and  specific  heat  for  seawater, 
respectively.  The  overbars  indicate  time-averaged  atmo¬ 
spheric  quantities. 

The  atmospheric  component  provides  a  total  of  six  fields 
to  the  ocean  component,  including  the  sea  level  pressure 
(SLP),  the  surface  wind  stress  in  the  x  and  y  directions 
(Eqs.  4  and  5),  the  total  heat  and  moisture  fluxes  (Eqs.  6 
and  7),  and  the  net  solar  radiation.  The  net  solar  radiation  is 
used  to  compute  the  diabatic  heating  contribution  in  the 
potential  temperature  equation  (not  shown).  The  SLP  is  also 
needed  as  input  for  the  pressure  calculations  in  the 
momentum  equations  of  the  ocean  model.  The  coupling 
frequency  between  atmosphere  and  ocean  is  user- 
determined,  typically  on  the  order  of  30  min  or  less. 

2.3  NAVDAS 

The  NRL  Atmospheric  Variational  Data  Assimilation  System 
(NAVDAS)  is  used  for  atmospheric  data  assimilation  in  the 
coupled  system.  NAVDAS  is  an  observation  space-based 
3DVAR  suite  for  generating  maximum  likelihood  atmospheric 
state  estimates  to  satisfy  a  variety  of  Navy  needs  ranging  from 
global  initial  conditions  for  Navy  global  prediction  to  local 
initial  conditions  for  forward-deployed  mesoscale  prediction 
(Daley  and  Barker  2001a).  It  is  used  by  Fleet  Numerical 
Meteorology  and  Oceanography  Center  (FNMOC)  in  the 
Navy  Operational  Global  Atmospheric  Prediction  System 
(NOGAPS)  and  the  operational  mesoscale  system 
(COAMPS)  as  their  data  assimilation  schemes. 

The  observation  space  form  of  the  3D  VAR  equation  is 
given  as  (Daley  and  Barker  2001a): 

X,-Xb=  PbHT(HPbHT  +  R T1  [y  -  ¥(*„)],  (8) 

where  xa  is  the  analysis,  xb  the  background  (forecast  or 
prior),  Pb  the  positive-definite  background  error  covariance 
matrix,  H  the  Jacobian  matrix  corresponding  to  the 
(possibly)  nonlinear  forward  or  observation  operator  'P,  R 
the  observation  error  covariance  matrix,  and  y  the  observa¬ 
tion  vector.  Here  xa-xb  is  the  correction  vector  (analysis 
increment)  and  y  —  ^(xb)  is  the  innovation  vector  (obser¬ 
vation  increment)  of  length  L.  The  matrix  to  be  inverted, 
HPbHT+R ,  is  an  L  x  L  symmetric  matrix.  A  complete 
description  of  the  formulation  and  development  of  NAVDAS 
can  be  found  in  Daley  and  Barker  (2001a). 


Springer 


1940 


Ocean  Dynamics  (2011)  61:1937-1954 


NAVDAS  directly  assimilates  a  variety  of  observation 
types  and  off-time  observations,  including  radiosondes  and 
pibals,  surface  observations  from  land  and  sea,  TIROS 
Operational  Vertical  Sounder  (TO VS)  radiances  or  tempera¬ 
ture  profiles,  Aircraft  Meteorological  Data  Reporting 
(AMDAR)  observations  and  pilot  reports,  and  Special  Sensor 
Microwave  Imager  (SSM/I)  surface  wind  speed  and  total 
precipitable  water.  Its  formulation  permits  considerable  local 
vertical  and  horizontal  anisotropy.  NAVDAS  features  vertical 
variation  of  the  horizontal  correlation  scales,  horizontal 
variation  of  the  vertical  correlation  lengths,  and  vertical 
variation  of  the  mass-divergent  wind  coupling. 

NAVDAS  interfaces  with  NOGAPS  (global)  and 
CO  AMPS  (regional)  through  the  correction  vectors  that  are 
produced  for  the  grid  of  that  model.  For  the  multiple  nested- 
grid  CO  AMPS,  the  same  set  of  observations  is  used  for  all 
nests  but  the  characteristics  of  the  background  error  covari¬ 
ance  vary  between  the  outer  grid  and  the  innermost  nest.  The 
observation  innovation  is  computed  from  the  highest  resolu¬ 
tion  background  field  and  the  analysis  grid  is  constmcted  by 
removing  any  duplicate  grid  points  from  the  CO  AMPS  grid 
meshes.  An  additional  correction  is  added  to  the  outer  grid 
meshes  to  account  for  differences  between  the  background 
forecasts  in  the  grid  meshes  (Sashegyi  et  al.  2009). 

One  of  the  advantages  in  using  NAVDAS  in  the  global  or 
regional  ensemble  forecast  system  is  that  it  produces  estimates 
of  the  error  variance  of  its  analyses.  Such  variance  estimates 
can  be  used  to  constrain  the  magnitude  of  initial  perturbations 
that  represent  ET  transformations  or  linear  combinations  of 
ensemble  forecast  perturbation  (Bishop  and  Toth  1999).  It 
overcomes  the  drawback  of  the  multivariate  optimum 
interpolation  (MVOI)  scheme  used  in  earlier  tests  of  the 
CO  AMPS  ensemble  prediction  system  since  the  MVOI 
scheme  does  not  produce  estimates  of  the  analysis  error 
variance  (Bishop  et  al.  2009;  Holt  et  al.  2009). 

NAVDAS  analysis  error  variance  estimates  are  used  in  the 
ET  scheme  for  the  NRL  global  ensemble  prediction  system 
and  results  show  that  the  generated  ET  analysis  perturbations 
exhibit  statistically  significant,  realistic  multivariate  correla¬ 
tions.  The  forecast  ensembles  are  comparable  to  or  better  in  a 
variety  of  measures  than  those  produced  by  FNMOC  bred- 
growing  modes  scheme  (McLay  et  al.  2008).  Even  the 
NAVDAS  analysis  error  variance  generated  from  the  global 
model  and  interpolated  for  the  mesoscale  ensemble  predic¬ 
tion  system  is  able  to  produce  reasonable  ET  analysis 
perturbations  (Bishop  et  al.  2009;  Holt  et  al.  2009). 

2.4  NCODA 

The  ocean  data  assimilation  component  of  the  coupled 
ensemble  forecast  system  is  the  Navy  Coupled  Ocean  Data 
Assimilation  (NCODA)  system.  NCODA  is  a  fully  multi¬ 
variate  3DVAR  system  that  provides  simultaneous  analyses 


of  five  ocean  variables:  temperature,  salinity,  geopotential, 
and  u-  and  v-vector  velocity  components.  The  horizontal 
correlations  are  multivariate  in  geopotential  and  velocity, 
thereby  permitting  adjustments  to  the  mass  fields  to  be 
correlated  with  adjustments  to  the  flow  fields.  The  velocity 
adjustments  (or  increments)  are  in  geostrophic  balance  with 
the  geopotential  increments,  which  in  turn  are  in  hydrostatic 
agreement  with  the  temperature  and  salinity  increments. 
NCODA  corrects  ocean  model  initial  conditions  using  a 
sequential  incremental  update  cycle.  It  has  been  cycled  with 
a  variety  of  ocean  forecast  models,  including  HYCOM 
(HYbrid  Coordinate  Ocean  Model)  and  NCOM,  which  are 
in  operational  use  at  the  Navy  forecasting  centers.  The 
system  can  be  run  global  or  regional,  where  it  supports  re¬ 
locatable,  multi-scale  analyses  on  nested  higher-resolution 
grids  using  a  3:1  nested  grid  ratio.  This  nesting  strategy  is 
of  particular  importance  in  Navy  applications  where  very 
high  resolution  is  required  in  a  rapid  environmental 
assessment  mode  of  operation. 

NCODA  assimilates  a  wide  variety  of  ocean  observation 
data  types  including  satellite  SST  retrievals  (AATSR, 
AMSR-E,  METOP-AVHRR,  GOES,  MSG,  MTSAT, 
NOAA-AVHRR),  sea  ice  concentration  retrievals  (SSM/I, 
SSMIS),  fixed,  drifting  buoy,  and  ship  in  situ  SSTs; 
temperature  and  salinity  profiles  (Argo,  CTD  casts,  fixed 
and  drifting  buoys,  gliders,  expendable  bathythermo¬ 
graphs);  satellite  altimeter  sea  surface  height  anomalies 
(ENVISAT,  GFO,  Jason,  Topex);  significant  wave  height 
retrievals  from  altimeters  and  wave  buoys;  and  velocity 
observations  (HF  radar,  surface  drifters,  acoustic  Doppler 
current  profilers,  gliders,  and  Argo  trajectories).  NCODA  is 
tightly  coupled  to  a  fully  automated  ocean  data  quality 
control  system  that  is  executed  in  real  time  at  the  Navy 
centers  (Cummings  2011). 

NCODA  is  an  oceanographic  implementation  of  the 
NAVDAS  algorithm  and  solves  the  observation  space  form 
of  the  3D  VAR  in  the  ocean  using  Eq.  8.  A  specification  of 
the  background  and  observation  error  covariances  in  the 
analysis  is  very  important.  Background  error  covariances 
are  separated  into  an  error  variance  and  a  correlation.  The 
correlations  are  separated  further  into  horizontal  and 
vertical  components.  Horizontal  correlation  length  scales 
vary  with  location  and  are  set  proportional  to  the  first 
baroclinic  Rossby  radius  of  deformation  (Chelton  et  al. 
1998),  which  varies  from  ~10  km  at  the  poles  to  ~200  km 
at  the  equator.  Flow  dependence  is  introduced  by  modifying 
horizontal  correlations  with  a  tensor  that  takes  into  account 
the  geopotential  height  separation  between  observations 
and  grid  locations.  The  flow  dependence  spreads  model 
data  differences  along,  rather  than  across,  geopotential 
contours  such  as  ocean  fronts  and  eddies.  Vertical  correla¬ 
tion  length  scales  are  computed  by  scaling  a  specified 
change  in  density  stability  criterion  with  background 
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vertical  density  gradients  according  to: 

K=pJ{.dp/dz)  (9) 

where  hY  is  the  vertical  correlation  length  scale,  ps  is  the 
specified  change  in  density  criterion  (-0.12  kg  m-3),  and 
dp/dz  is  the  local  vertical  density  gradient.  The  resulting 
vertical  correlation  length  scales  vary  with  location  and 
depth  and  are  large  (small)  when  the  water  column 
stratification  is  weak  (strong).  Vertical  correlations  are 
updated  based  on  the  forecast  valid  at  the  analysis  time. 
In  this  way,  the  vertical  scales  evolve  from  one  analysis 
cycle  to  the  next  and  capture  changes  in  mixed  layer  and 
thermocline  depths. 

Background  error  variances  are  derived  from  a  time 
history  of  the  analysis  increment  fields  and  a  model 
variability  component  derived  from  a  time  history  of  model 
forecast  differences  according  to: 

m  m 

eb  =  V  Wj(xj  -  xj-i)2  +  ^r-aj  (10) 

j=  1  7=1 

where  eh  is  the  background  error  variance,  j  is  the  update 
cycle,  m  is  the  number  of  update  cycles  into  the  past,  Wj  is  a 
weight  computed  from  a  geometric  time  series,  Wj=(  l-</>y_1, 
where  0  is  a  tunable  constant  between  0  and  1  (typically  set 
to  0.1),  (Xj—Xj-i)  is  the  difference  between  successive 
forecasts,  aj  is  the  analysis  increment  field,  and  r  is  a 
correlation  time  scale  computed  from  a  ratio  of  the  spatially 
varying  horizontal  length  scales  h  and  forecast  velocity  fields 
v  valid  at  the  analysis  time,  i=e~ih/v\  The  weights  are 
normalized  such  that  the  weighted  averages  of  the  two 
components  are  unbiased.  Subscripts  depicting  3D  position 
coordinates  have  been  eliminated  for  clarity.  The  time  scales 
vary  with  location  and  depth  and  range  from  days  near  the 
surface  in  the  western  boundary  current  regions  to  years  at 
depth  in  the  tropical  ocean  where  current  speeds  are  slow 
and  length  scales  are  large.  In  this  approach,  model  data 
errors  tend  to  dominate  the  background  error  variance 
estimates  in  low-flow  environments,  while  model  variability 
tends  to  dominate  in  high-flow  (i.e.,  hurricane-induced) 
environments.  Background  error  variances  are  computed 
for  all  of  the  model  prognostic  variables  and  updated  at  each 
analysis  update  cycle. 


3  Description  of  the  Ensemble  Transform  technique 

The  ET  ensemble  generation  technique  transforms  an 
ensemble  of  forecast  perturbations  into  an  ensemble  of 
analysis  perturbations  that  are  consistent  with  a  user- 
provided  estimate  of  the  analysis  error  covariance  matrix 
Pa.  Assume  that  there  are  n  model  variables  in  total 
(typically,  n  ~  O(107))  and  K  ensemble  members  in  total 


(typically,  K  ~  0(10)  —  O(102)).  Let  the  columns  of  the  n 
x  K  matrix  list  the  forecast  perturbations  about  the 
forecast  ensemble  mean  all  valid  at  the  time  of  the  most 
recent  analysis.  Let  the  columns  of  the  n  x  K  matrix  Xa 
denote  the  analysis  perturbations  that  will  be  added  to  the 
most  recent  analysis  to  create  an  ensemble  of  initial 
conditions  consistent  with  the  accuracy  of  the  most  recent 
analysis.  The  ET  transforms  Xf  into  Xa  using: 

xa  =  xy  (ii) 


As  noted  in  McLay  et  al.  (2008),  to  compute  the  required 
transformation  matrix  one  first  performs  the  eigenvector 
decomposition: 


X/rP“-1X/ 


n 


=  CMC7-, 


(12) 


where  the  orthonormal  K  x  K  matrix  C  lists  the  orthonor¬ 
mal  eigenvectors  of  X,T  Pa_1  xfln  and  the  diagonal  matrix 
A  lists  the  corresponding  eigenvalues.  Because  the  sum  of 
the  forecasts  perturbations  is  zero,  one  of  the  eigenvalues  will 
always  be  equal  to  zero.  The  eigenvector  corresponding  to  this 
eigenvalue  is  always  parallel  to  the  vector  ^-elements 


ir  =  i] 


where  1  is  a  column  vector  of  ones.  For  the  reasons 
discussed  in  McLay  et  al.  (2008),  a  new  diagonal  matrix 
A  is  created  from  A  by  replacing  the  zero  eigenvalue  by  the 
unity.  The  transformation  matrix  is  then  given  by: 


T  =  C4-1/2Cr, 


(13) 


When  T  is  defined  in  this  way,  the  analysis  perturbations 
obtained  from  Eq.  1 1  satisfy  the  equation: 


XaTPa-lXa  _  llr 
n 


(14) 


where  i  is  a  K  x  K  identity  matrix.  As  noted  in  McLay  et  al. 
(2008),  Eq.  14  implies  that  (1)  the  analysis  perturbations  are 
quasi-orthogonal  under  the  analysis  error  covariance  norm 
and  (2)  if  the  ensemble  size  K  was  large  enough  to  satisfy 
K=n,  then  the  sample  covariance  of  the  analysis  perturba¬ 
tions  would  equal  the  prescribed  analysis  error  covariance 
matrix  Pa.  Note  that  Eq.  11  means  that  the  analysis 
perturbations  are  just  linear  combinations  of  forecast  pertur¬ 
bations.  Hence,  if  the  forecast  perturbations  lie  in  the  space 
of  growing  balanced  perturbations,  then  the  analysis  pertur¬ 
bations  will  also.  Indeed  Eq.  1 1  means  that  the  ET  ensemble 
generation  technique  may  be  viewed  as  an  advanced  form  of 
Toth  and  Kalnay’s  (1993)  breeding  technique  with  added 
constraints  of  quasi-orthogonality  and  fit  to  the  prescribed 
analysis  error  covariance  matrix  Pa. 

In  its  purest  form,  only  one  set  of  transformation 
coefficients  is  used  for  the  entire  state.  However,  McLay 
et  al.  (2010)  have  shown  that  smoothly  blending  the 
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transformation  coefficients  relevant  to  large  regions  such  as 
the  tropics  and  extra-tropics  can  lead  to  performance 
advantages.  In  this  first  study,  we  approximate  Eq.  12  for 
both  the  atmospheric  and  oceanic  components  and  apply 
separate  transformation  coefficients  for  the  atmosphere  and 
the  ocean. 

For  the  ocean  component,  NCODA  includes  options  to 
estimate  the  analysis  error  variance  defined  as  the  reduction  of 
forecast  error  from  assimilation  of  the  observations.  Similar  to 
the  description  of  NAVDAS  above  (“Section  2.3”)  the 
primary  application  of  the  NCODA  analysis  error  variance 
is  in  the  oceanic  ET  technique.  The  ocean  ET  is  multivar¬ 
iable  and  is  computed  for  temperature,  salinity,  and  velocity 
simultaneously.  In  warm  start  mode,  the  ocean  ET  operates 
identical  to  the  atmosphere  ET  (see  “Section  4”  for  a 
discussion  of  cold  start  and  warm  start  modes  of  the  coupled 
system).  However,  in  cold  start  mode,  the  ocean  ET  is 
limited  by  two  factors:  (1)  there  is  no  forecast  error  product 
from  the  deterministic  global  NCOM  (gNCOM)  and  (2) 
there  is  no  global  NCOM  ensemble  from  which  to  draw  the 
initial  set  of  perturbations.  To  solve  this  problem  in  cold  start 
mode,  we  treat  a  time  series  of  gNCOM  forecasts  at  the 
update  cycle  interval  as  members  of  an  ensemble.  The 
number  of  members  used  in  the  actual  ensemble  determines 
how  far  back  in  time  we  have  to  obtain  gNCOM  fields.  The 
NCODA  analysis  background  error  variances  and  the  initial 
perturbations  for  the  ensemble  are  initialized  from  this  time 
series.  This  process  results  in  NCODA  analysis  background 
error  variances  and  ET  perturbations  that  are  multivariate, 
balanced,  and  consistent.  The  ability  to  relocate  the  coupled 


system  for  rapid  environmental  assessment  is  thereby 
enhanced  as  well. 


4  Implementation  of  the  mesoscale  ensemble  system 

The  components  of  the  ocean-atmosphere  coupled  system 
as  described  in  “Section  3”  are  efficiently  integrated  into  an 
ensemble  framework  under  the  Earth  System  Modeling 
Framework  as  illustrated  in  Fig.  1.  A  basic  differential 
starting  point  for  the  ensemble  system  is  a  “cold  start” 
versus  a  “warm  start.”  The  specifics  of  each,  and  their 
differences,  are  described  here. 

A  cold  start  is  executed  initially  and  only  once  for  a 
given  case  study.  The  first  module  executed  in  the  system  is 
the  analysis,  indicated  by  the  upper  box  in  Fig.  la.  Because 
this  is  the  initial  execution  of  the  mesoscale  system,  there 
are  no  previous  “background”  (BKG)  forecasts  to  be  used 
as  first  guesses  for  the  analysis.  Thus,  global  fields  must  be 
used.  The  control  member  is  the  only  ensemble  member  to 
execute  the  analysis  (step  1  in  Fig.  1).  The  analysis  consists 
of  two  basic  parts:  preparing  the  observations,  initial 
conditions  (IC),  boundary  conditions  (BC),  and  BKG  fields 
for  (1)  the  atmosphere  and  (2)  for  the  ocean.  The 
atmospheric  analysis  is  based  on  NAVDAS  (see  “Sec¬ 
tion  2.3”),  using  the  control  member  from  the  3 3 -member 
(32  members  plus  control)  Navy  Operational  Global 
Atmospheric  Prediction  System  (NOGAPS)  ensemble  run 
at  T119  horizontal  resolution  (~110  km)  and  30  vertical 
levels,  interpolated  to  the  COAMPS  atmospheric  domain  as 
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Fig.  1  Schematic  showing  the  implementation  of  the  air-ocean  coupled  ensemble  system  for  a  a  cold  start  and  b  a  warm  start 
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the  first  guess.  The  oceanic  analysis  is  based  on  NCODA 
(see  “Section  2.4”),  using  gNCOM  interpolated  to  the 
COAMPS  oceanic  domain  as  first  guesses.  There  is 
currently  no  global  ocean  ensemble,  so  the  deterministic 
run  is  used. 

The  analyzed  atmospheric  and  oceanic  output  fields 
from  the  control  member  are  next  copied  to  all  members 
(step  2).  Each  member  now  creates  unique  atmospheric 
lateral  BC  from  28  members  of  the  global  NOGAPS 
ensemble  (step  3).  The  members  of  the  global  ensemble 
are  produced  using  perturbation  methodology  based  on  the 
banded  ET  as  described  in  McLay  et  al.  (2008,  2010).  The 
banded  ET  produces  initial  perturbations  by  transforming  6- 
h  ensemble  forecast  perturbations  to  be  consistent  with 
analysis  error  estimates.  The  local  banded  ET  performs  this 
transformation  in  latitude  bands,  resulting  in  a  better  match 
to  the  analysis  error  estimate,  and  much  improved  perfor¬ 
mance  in  the  mid-latitudes,  as  compared  to  the  global  ET 


computed  over  the  entire  globe.  Stochastic  kinetic  energy 
backscatter  (Shutts  and  Palmer  2004;  Shutts  2005;  Berner 
et  al.  2009)  is  planned  for  future  implementations  of  the  ET 
adopted  to  account  for  model  uncertainty  and  increase  in 
size  of  tropical  perturbations.  Step  4  is  to  execute  the 
coupled  forecast  model  for  each  of  the  members,  using  the 
same  initial  conditions  but  different  lateral  atmospheric  BC, 
followed  by  post-processing  (step  5). 

The  implementation  of  the  ensemble  system  for  a  warm 
start  differs  from  a  cold  start  in  several  ways  (Fig.  lb).  First, 
though  the  analysis  is  again  executed  only  for  the  control 
member  (upper  box  in  Fig.  lb),  the  BKG  fields  for  both 
NAVDAS  and  NCODA  are  previous  control  member 
forecasts  (typically  6  or  12  h)  from  the  ensemble  (not  global 
fields).  The  analysis  creates  and  copies  the  same  fields  as  for 
the  cold  start  (step  2).  However,  the  IC  for  the  individual 
members  are  now  obtained  by  adding  the  ET  perturbations  to 
the  control  analysis  (see  “Section  3”),  separate  for  both  the 


Table  1  Listing  of  the  atmospheric  physics  options  and  perturbations.  Italicized  values  for  the  members  (1-28)  indicate  values  that  differ  from  the 
control  (member  0;  see  text  for  a  detailed  description  of  the  physics  parameters) 


Member 

abl 

mixlen 

Flux 

w-kf 

tinc-lcl 

eld-rad 

precip 

Graupel 

Auto-conv 

Rain-int 

Snow-int 

0  (control) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

1  (PBL  buoyancy) 

1 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

2  (PBL  mixing  length) 

2 

1.25 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

3  (PBL  mixing  length) 

2 

0.75 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

4  (PBL  mixing  length) 

2 

1.5 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

5  (PBL  mixing  length) 

2 

0.5 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

6  (PBL  surface  flux) 

2 

1.0 

1.25 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

7  (PBL  surface  flux) 

2 

1.0 

0.75 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

8  (PBL  surface  flux) 

2 

1.0 

1.5 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

9  (PBL  surface  flux) 

2 

1.0 

0.5 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

10  (cumulus  W  trigger) 

2 

1.0 

1.0 

1.5 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

11  (cumulus  W  trigger) 

2 

1.0 

1.0 

0.5 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

12  (feedback  fraction) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.5 

True 

0.0004 

8.0e6 

2.0e7 

13  (cumulus  T  trigger) 

2 

1.0 

1.0 

1.0 

1.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

14  (cumulus  T  trigger) 

2 

1.0 

1.0 

1.0 

-1.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

15  (cumulus  T  trigger) 

2 

1.0 

1.0 

1.0 

2.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

16  (cumulus  T  trigger) 

2 

1.0 

1.0 

1.0 

-2.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

17  (cloud  radius) 

2 

1.0 

1.0 

1.0 

0.0 

500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

18  (cloud  radius) 

2 

1.0 

1.0 

1.0 

0.0 

1,000.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

19  (cloud  radius) 

2 

1.0 

1.0 

1.0 

0.0 

2,000.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

20  (cloud  radius) 

2 

1.0 

1.0 

1.0 

0.0 

3,000.0 

0.0 

True 

0.0004 

8.0e6 

2.0e7 

21  (feedback  fraction) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

1.0 

True 

0.0004 

8.0e6 

2.0e7 

22  (graupel) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

False 

0.0004 

8.0e6 

2.0e7 

23  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.00004 

8.0e6 

2.0e7 

24  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.004 

8.0e6 

2.0e7 

25  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8. 0e7 

2.0e7 

26  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8. 0e5 

2.0e7 

27  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0  e6 

2.0e8 

28  (microphysics) 

2 

1.0 

1.0 

1.0 

0.0 

1,500.0 

0.0 

True 

0.0004 

8.0e6 

2.0e6 
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Table  2  Description  of  the  TC  case  studies.  The  number  of  forecasts  indicates  the  number  of  data  assimilation  forecasts  executed  every  12  h  for 
the  case  study.  Ike-u  is  the  same  as  Ike  except  that  the  atmospheric  model  was  run  uncoupled  with  the  ocean  model 


Case 

Cold  start 
date 

No.  of 
forecasts 

Members 

Atmos  resolution 
(km) 

Ocean  resolution 
(km) 

Atmospheric 

levels 

Ocean 

levels 

Forecast 
length  (h) 

Moving 

nests 

Coupled 

Hanna 

2008090112 

6 

29 

81,  27 

27 

30 

40 

54 

Y 

Y 

Ike 

2008090500 

14 

29 

81,  27 

27 

30 

40 

60 

Y 

Y 

Ike-u 

2008090500 

14 

29 

81,  27 

27 

30 

40 

60 

Y 

N 

atmosphere  and  the  ocean  (step  3).  The  lateral  boundary 
conditions  are  created  the  same  as  for  a  cold  start  (step  4)  as 
are  the  execution  of  the  forecast  model  (step  5)  and  post¬ 
processing  (step  6).  The  entire  process  (from  step  1  to  6)  is 
repeated  for  the  next  data  assimilation  cycle. 

Previous  simulations  testing  the  ensemble  system 
displayed  tendencies  for  the  forecast  to  be  under- 
dispersive  (too  little  spread).  To  address  this  issue  and 
to  better  represent  uncertainty,  atmospheric  perturbations 
are  developed  and  tested  as  described  in  the  following 
section. 

4.1  Atmospheric  perturbations 

In  order  to  represent  the  uncertainty  due  to  the  model 
physical  parameterizations,  each  of  the  CO  AMPS  mem¬ 
bers  is  run  with  a  different  set  of  physical  parameteriza¬ 
tion  values  as  shown  in  Table  1.  Previous  COAMPS 
ensemble  experiments  were  conducted  with  a  variety  of 
atmospheric  parameters.  The  results  from  these  experi¬ 
ments  isolated  several  key  parameters  in  the  ABL  and 
cumulus  parameterizations  and  cloud  microphysics  to 
which  COAMPS  forecasts  were  most  sensitive.  The 
parameterization  developers  were  consulted  in  order  to 
focus  on  a  subset  of  key  parameters  that  they  thought 
represented  the  largest  uncertainty.  These  parameters  are: 
(1)  the  type  of  ABL  parameterization,  (2)  the  magnitude  of 
the  vertical  mixing  length,  (3)  the  magnitude  of  the  surface 
sensible  and  latent  heat  fluxes,  (4)  the  grid-scale  vertical 
velocity  used  in  the  K-F  cumulus  parameterization  trigger 
for  convective  initiation,  (5)  the  fraction  of  available 
convective  precipitation  partitioned  to  the  grid-scale 
precipitation,  (6)  the  temperature  increment  at  the  lifted 
condensation  level  (LCL)  (which  also  impacts  the  con¬ 
vective  initiation  in  the  K-F  parameterization),  (7)  the 
cloud  updraft  radius  used  in  the  K-F  parameterization,  (8) 
microphysical  interactions  that  include  graupel,  (9)  the 
microphysics  auto-conversion  factor,  and  (10)  the  micro¬ 
physics  slope  intercept  factors  for  snow  and  rain.  A  brief 
description  for  each  is  given  here. 

1.  The  ABL  parameterization:  Ensemble  member  1  uses 
the  original  Mellor  and  Yamada  (1974,  1982)  ABL 


scheme  (abl=l).  The  modified  Mellor  and  Yamada 
boundary  layer  scheme  (abl=2)  has  several  features  not 


Fig.  2  The  coupled  model  domain  setup  for  the  TC  case  studies: 
a  Hanna  and  b  Ike  (see  Table  2  for  details  on  each  configuration) 
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available  in  the  original  formulation.  A  statistical  cloud 
fraction  is  computed  for  the  buoyancy  calculation.  This 
is  used  in  diagnosing  the  virtual  buoyancy  flux  which 
impacts  the  buoyant  production  of  turbulence.  The 
thermodynamic  quantities  used  for  determining  stabil¬ 
ity  include  liquid  water  potential  temperature  and  a 
hybrid  potential  temperature  incorporating  both  liquid 
water  potential  temperature  and  virtual  potential  tem¬ 
perature.  All  five  species  of  water  substance  (vapor, 
liquid,  snow,  ice,  and  graupel)  are  used  in  all 
thermodynamic  computations.  The  use  of  these  quan¬ 
tities  removes  the  discontinuity  at  cloud  base  which 
was  a  common  occurrence  in  the  original  Mellor  and 
Yamada  scheme.  In  terms  of  shear  production,  a 
background  shear  based  on  the  turbulence  level  is 
included  as  the  lower  limit  on  shear  production  of  TKE. 
The  turbulence  length  scale  is  modified  to  increase 
more  rapidly  with  elevation  in  highly  convective 
situations.  In  addition,  the  Brunt- Vaisala  length  scale 
used  in  stable  conditions  is  not  utilized  within  cloud 
layers.  The  modified  ABL  scheme  is  used  for  all 
members  except  ensemble  member  1  (see  Table  1). 


2.  The  magnitude  of  the  vertical  mixing  length:  CO  AMPS 
uses  separate  vertical  eddy  mixing  coefficients  for 
momentum,  KM : 

=  SWv£-1^25  (15) 

and  for  the  scalar  variables,  Kn: 

Kli=Slihe-1'2,  (16) 

where  /v  represents  the  vertical  mixing  length,  SM  and 
Su  the  stability  factors  for  momentum  and  scalars,  and 
e  the  TKE  (Mellor  and  Yamada  1974;  Therry  and 
Lacarrere  1983).  This  vertical  mixing  length  is  scaled 
by  factor  mixlen  (ranging  from  0.5  to  1.5)  for 
perturbation  members  2-5. 

3.  The  magnitude  of  the  surface  sensible  and  latent  heat 
fluxes:  The  sensible  and  latent  fluxes  as  given  by 
Eqs.  2  and  3  are  scaled  by  factor  flux  (ranging  from 


Surface  v-current  (m  s1) 


Hanna:20080901 1 2 


10-m  v-wind  (m  s'1) 


New  interpolated  moving  nest  region  for 
atmospheric  ET 


5 


0 


L±5 


Fig.  3  The  ocean  nest  1  (27  km)  analyzed  surface  v-current  (m  s_1) 
and  the  location  of  the  atmospheric  nest  2  (27  km)  for  the  control 
member  ( solid  box)  and  member  1  (moved  nest  indicated  by  dashed 


lines)  for  the  1200  UTC  1  September  Hurricane  Hanna  case  (left)  and 
the  analyzed  nest  2  (27  km)  10-m  wind  v-component  (m  s-1)  after  the 
moving  nest  ET  algorithm  has  been  applied  (right) 
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Hanna:  2008090112-2008090700 


Ocean  prediction  velocity  errors  (cm  s_1) 


0  5  10  15  20  25  30 


Fig.  4  Ocean  prediction  velocity  errors  (cm  s  l)  for  ocean  nest  1  (27-km)  for  Hurricane  Hanna  averaged  from  1200  UTC  1  September  to  0000 
UTC  7  September  for  depths  a  1  m,  b  20  m,  c  40  m,  d  85  m,  e  130  m,  and  f  200  m 


0.5  to  1.5)  over  land  and  water  for  perturbation 
members  6-9. 

4.  The  grid-scale  vertical  velocity  used  in  the  K-F 
cumulus  parameterization  to  diagnose  convective  initi¬ 
ation:  The  value  is  scaled  by  factor  w-kf  (ranging  from 
0.5  to  1.5)  for  members  10  and  11. 

5.  The  fraction  of  available  convective  precipitation  that  is 
partitioned  to  the  grid- scale  precipitation:  The  value  is 
scaled  by  factor  precip  (ranging  from  0.5  to  1.0)  for 
members  12  and  21. 

6.  The  temperature  increment  at  the  LCL:  The  increment  ( K ) 
added  to  the  temperature  at  the  LCL  (i tinc-lcl )  ranges 
from  -2.0  to  2.0  K  for  members  13-16.  This  increment 
has  an  impact  on  the  convective  initiation. 

7.  The  cloud  updraft  radius  used  in  the  K-F  parameteri¬ 
zation:  The  radius  cld-rad  (m)  varies  from  500  to 
3,000  m  for  members  17-20. 

8.  Graupel:  Graupel  microphysics  calculations  ( graupel ) 
are  not  included  for  member  22. 

9 .  Microphysics  auto-conversion  factor:  The  cloud  water  mixing 
ratio  threshold  in  the  Kessler  (1969)  auto-conversion  of  cloud 
water  to  rain  (< auto-conv )  varies  from  0.00004  to  0.004  for 


Hanna  Nest  2:  TC  forecast  track  from  2008090400 


Fig.  5  TC  forecast  tracks  for  29-member  Hurricane  Hanna  for  the 
warm  start  at  0000  UTC  4  September  for  0  to  54  h  for  the  COAMPS 
coupled  ensemble.  The  white  circles  on  the  best  track  indicate  the  day 
(04,  05,  06),  and  the  solid  circles  for  the  forecast  tracks  are  every  12  h. 
(Note  that  the  member  number  colors  are  not  consistent  for  the  map 
and  the  error  statistics  are  given  in  Fig.  6) 
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Fig.  6  COAMPS  coupled  en¬ 
semble  statistics  for  29-member 
Hurricane  Hanna  for  the  warm 
start  at  0000  UTC  4  September 
for  0  to  54  h  of  a  maximum 
wind  speed  (kts),  b  track  error 
(nm),  and  c  maximum  wind 
speed  error  (kts).  The  control 
forecast  is  given  by  the  thick 
black  line ,  and  the  ensemble 
mean  is  given  by  the  thick  blue 
line.  The  wind  speed  from  the 
observed  best  track  TC  is  given 
by  the  solid  line.  (Note  that  the 
member  number  colors  are  not 
consistent  for  the  map  in  Fig.  5 
and  the  error  statistics) 
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members  23  and  24.  There  is  no  conversion  unless  cloud 
water  mixing  ratio  >  auto-conv. 

10.  Microphysics  slope  intercept  factors:  The  intercept 
factors  for  the  raindrop  size  (rain-int)  vary  from  8.e5 
to  8.e7  for  members  25  and  26  and  snow  flake  size 
( snow-int )  distributions  vary  from  2.e6  to  2.e8  for 
members  27  and  28. 


5  Results  of  tropical  cyclone  case  studies 

Coupled  ensemble  simulations  are  conducted  for  two  TC 
case  studies:  Hanna  (1-6  September  2008)  and  Ike  (5-14 
September  2008).  The  primary  focus  of  this  study  is  to 
demonstrate  the  capability  of  the  coupled  ensemble  system 
not  to  provide  high-resolution  simulations  of  the  TCs. 
Thus,  the  following  analyses  stress  the  performance  of  the 
coupled  ensemble  (particularly  the  ensemble  mean  versus 
the  deterministic  control  forecast)  and  not  the  detailed  high- 
resolution  validation  against  observations. 

Hurricane  Hanna  is  investigated  over  the  period  from  1 
to  6  September  2008.  However,  the  storm  was  prominent  as 
early  as  28  August  when  it  was  named.  By  30  August,  the 
circulation  was  weak,  centered  near  its  location  just  east  of 
the  Bahamas.  Hurricane  Gustav’s  circulation  influenced 
Hanna’s  track  on  1  September,  causing  a  slow  drift  to  the 
south-southwest  with  convection  and  strong  intensification, 
with  the  National  Hurricane  Center  (NHC)  upgrading 
Hanna  to  a  hurricane  (sustained  winds  of  at  least  64  kts  or 
33  m  s_1).  However,  Gustav  again  impacted  Hanna’s 
structure  and  by  2  September  strong  wind  shear  weakened 
the  circulation  and  the  storm  was  downgraded  to  a  tropical 
storm  on  3  September.  Hanna  made  landfall  near  the  South 
Carolina-North  Carolina  border  on  6  September  and 
subsequently  went  extra-tropical  by  7  September. 

Hurricane  Ike  is  investigated  over  the  period  from  5  to 
14  September  2008.  Ike  was  a  very  long-lived  system,  first 
evident  in  the  Atlantic  as  Tropical  Depression  Nine  on  1 
September,  with  several  subsequent  periods  of  intensifica¬ 
tion  and  weakening  over  the  next  10  days.  Ike  reached  peak 
intensity  (—145  kts  or  79  m  s_1)  by  4  September,  weakened 
to  category  3  on  5  September  and  category  2  on  6 
September.  This  weakening  was  short-lived  and  rapid 
intensification  occurred  again  in  the  early  hours  of  7 
September,  reaching  category  4  (with  sustained  winds  of 
115  kts  near  the  Caicos  Islands)  just  6  hours  after  being 
downgraded  to  category  2.  Ike  made  several  landfalls  at 
Caribbean  islands  on  8  and  9  September  before  rapidly 
intensifying  during  the  night  of  1 0  September  (reduction  in 
central  SLP  from  963  to  944  hPa  and  increase  in  wind 
speed  from  74  to  87  kts)  as  it  passed  over  the  loop  current 
in  the  Gulf  of  Mexico.  From  10  to  12  September,  Ike 


remained  a  category  2  storm  (-95  kts)  as  it  moved  north- 
northwestward  toward  Galveston,  TX,  USA,  but  with  an 
unusually  large  spatial  distribution  of  high  wind  speeds, 
resulting  in  a  large  storm  surge  as  it  made  landfall  at  -0700 
UTC  13  September. 

A  variety  of  configurations,  including  differing  model 
resolution,  number  of  members,  forecast  length,  and 
moving  and  fixed  nest  configurations,  have  been  tested 
for  the  coupled  ensemble  system.  A  “standard”  configura¬ 
tion  is  chosen  for  the  two  TC  cases  presented  here,  as  listed 
in  Table  2  and  shown  in  Fig.  2.  Each  of  the  studies  begins 
with  a  cold  start  (as  described  in  “Section  4”),  employs  data 
assimilation  every  12  h,  and  uses  atmospheric  model 


Fig.  7  TC  forecast  tracks  for  29-member  Hurricane  Ike  for  the  warm 
start  at  0000  UTC  11  September  for  the  COAMPS  a  coupled 
ensemble  and  b  uncoupled  ensemble.  The  tracks  correspond  to 
different  forecast  lengths  (60  and  72  h)  than  the  best  track  ( solid 
black).  The  white  circles  on  the  best  track  indicate  the  day  (11,  12,  13) 
and  the  solid  circles  for  the  forecast  tracks  are  every  12  h 
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Fig.  8  COAMPS  coupled  en¬ 
semble  statistics  for  29-member 
Hurricane  Ike  for  the  warm  start  at 
0000  UTC  1 1  September  for  0  to 
54  h,  similar  to  Fig.  6 
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perturbations  (as  outlined  in  “Section  4.1”),  with  a  two-way 
exchange  of  information  every  time  step  between  the  grid 
meshes.  The  Hanna  case  study  cold  start  begins  at  1200 
UTC  1  September  2008  (Ike  at  0000  UTC  5  September 
2008),  with  a  series  of  54-h  (60-h  for  Ike)  forecasts  every 
12  h  for  3  days  or  six  total  forecasts  (7  days  or  14  forecasts 
for  Ike).  The  ensemble  simulations  use  two  atmospheric 
domains  with  a  horizontal  resolution  of  8 1  and  27  km,  with 
the  second  nest  moving  automatically  with  the  TC  center, 
and  one  fixed  ocean  domain  with  27-km  resolution.  The 
atmospheric  model  has  30  vertical  levels  with  the  vertical 
grid  spacing  varied  from  10  m  near  the  surface  to  a 
maximum  of  7,500  m  near  the  model  top  (~30  km).  The 
ocean  model  has  40  vertical  levels,  including  19  terrain¬ 
following  sigma  levels  in  the  upper  140  m  and  21  fixed- 
depth  levels  between  140  m  and  the  ocean  bottom.  The  Ike 
case  study  is  also  simulated  with  the  atmospheric  model 
uncoupled  to  the  ocean  model  (Ike-u)  in  addition  to  the 
standard  coupled  configuration. 

For  TC  applications  in  which  the  inner  nests  of  the 
atmospheric  model  may  move  with  the  TC,  there  is  an 
inherent  problem  in  implementing  the  ET  to  create  a  new 
initial  state  because  it  cannot  be  guaranteed  that  the  model 
domains  of  every  member  will  coincide  spatially.  In  fact,  it 
is  a  desired  result  for  none  of  the  members  to  coincide,  such 
that  the  TC  track  of  each  member  is  different  to  enable  the 
spread  to  better  represent  the  natural  variability.  For  the 
application  of  the  ET,  this  results  in  a  spatially  inconsistent 
background  field.  Thus,  a  new  capability  to  apply  the  ET 
for  moving  nests  has  been  developed.  The  new  algorithm  is 
described  below. 

For  a  given  BKG  field  (i.e.,  12-h  forecast),  the  spatial 
area  for  each  member  that  is  common  with  the  control 
member  (i.e.,  overlap  region)  is  computed.  Then,  for  each 
individual  member,  that  area  of  the  moving  nest  which 
is  outside  the  overlap  region  is  interpolated  from  the 
domain  of  the  parent  (i.e.,  if  nest  2  is  the  moving  nest, 
then  nest  1  would  be  the  parent  used  for  interpolation). 
This  new  domain  that  is  now  the  same  as  the  control 
domain  is  used  as  the  BKG  field  for  the  ET.  Figure  3 
shows  an  example  of  the  implementation  of  the  algorithm 
for  the  1200  UTC  1  September  2008  Hanna  case.  The 
atmospheric  domain  for  the  control  member  is  given  by 
the  solid  box  in  the  right  corner  of  the  ocean  nest  1 
(27  km)  surface  v-current  (left  panel).  The  moving  nest 
member  has  moved  to  the  northwest  following  the  TC,  so 
the  new  interpolated  region  encompasses  a  region  to  the 
north  and  west  of  the  control  domain,  indicated  by  the 
dashed  line.  The  new  atmospheric  10-m  v-wind  compo¬ 
nent,  including  the  interpolated  region,  is  given  in  the 
right  panel.  Note  that  the  interpolated  region  contains 
only  weak  gradients.  The  amount  of  discontinuity  will 
obviously  depend  upon  the  difference  in  the  model 


solution  on  the  parent  and  child  domains  in  the 
interpolated  region.  In  addition,  the  smaller  the  overlap 
region  of  different  member  domains  (i.e.,  the  larger  the 
interpolated  regions),  the  larger  the  potential  adverse 
impact  on  the  ensemble  mean  because  of  increased 
interpolation. 

Figure  4  shows  the  NCODA  3D  VAR  ocean  velocity 
BKG  errors  for  the  Hurricane  Hanna  case  study  averaged 
from  1200  UTC  1  September  to  0000  UTC  7  September. 
As  discussed  in  “Section  2.4,”  there  are  two  components  of 
the  BKG  error:  model  variability  and  model  data  errors. 
Model  variability  is  dependent  on  the  time  history  of 
forecast  differences  at  the  updated  cycle  interval  (12  h),  and 
model  data  errors  are  determined  by  the  time  history  of  the 
analyzed  increment  fields.  As  a  result,  NCODA  3D  VAR 
BKG  errors  adapt  as  the  ocean  model  evolves  over  time  due 
to  changes  in  the  atmospheric  forcing  and  initial  conditions 
from  the  assimilation.  The  characteristics  of  the  ocean 
prediction  errors  of  temperature  and  salinity  (not  shown) 
have  similar  features  as  ocean  velocity.  For  hurricane- 
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Fig.  9  Mean  and  spread  track  statistics  for  29-member  Hurricane  Ike 
(coupled)  and  Ike-u  (uncoupled)  simulations  for  0000  UTC  11 
September.  The  track  error  values  are  plotted  to  54  h  to  correspond 
approximately  to  landfall  time 
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induced  high- flow  regimes,  Fig.  4  shows  that  the  horizontal 
distribution  of  the  ensemble  mean  ocean  velocity  BKG 
errors  within  the  ocean  mixed  layer  is  mainly  in  the  vicinity 
of  the  envelope  of  ensemble  tracks  near  the  hurricane 
maximum  winds. 

On  4  September,  Hanna  moved  out  of  a  counter¬ 
clockwise  looping  motion  over  the  Caicos  Island  and  began 
to  move  northwestward  (Fig.  5).  Hanna  turned  northward 
24  to  36  h  later  along  the  western  edge  of  a  now  better 
established  subtropical  ridge  over  the  Atlantic.  The  inten¬ 
sity  of  Hanna  remained  basically  unchanged,  around  55- 
60  kts  during  this  period  (Fig.  6a).  In  general,  the  ensemble 
mean  forecast  track  for  Hanna  follows  the  best  track  closely 
for  the  first  36  h  but,  over  the  next  18  h,  shifts  to  the  west 
of  the  best  track  as  the  storm  separates  from  an  upper-level 
low  pressure  system  and  is  influenced  by  the  westward 
extension  of  the  ridge.  A  similar  deviation  of  the  forecast 
track  is  also  evident  in  the  NOGAPS  forecasts  (not  shown). 
The  ensemble  mean  track  errors  (Fig.  6b)  at  24,  36,  and 
48  h  are  well  within  the  average  errors  for  the  official  NHC 
forecasts  of  Hanna  [1  nautical  mile  (nm)  =  1.852  km].  The 
ensemble  forecasts  show  average  intensity  errors  of  <10  kts 

Fig.  10  Sea  surface  temperature 
difference  (°C)  from  0000  UTC 
11  September  to  0000  UTC  13 
September  for  Hurricane  Ike  for 
a  TRMM  Microwave  Imager 
satellite  observations  and  b 
coupled  ensemble  mean.  The 
contours  in  b  are  sea  surface 
height  to  show  the  location  of 
the  warm  core  eddy  in  the  Gulf 
of  Mexico.  Only  surface  heights 
greater  than  0.4  m  are  drawn. 

The  ensemble  mean  track  is 
computed  from  the  0-  to 
11-h  average  forecast  position 


for  the  period  prior  to  landfall  (Fig.  6c),  while  official 
forecasts  over-estimate  Hanna  as  a  category-2  storm  upon 
landfall.  Furthermore,  the  ensemble  mean  of  intensities 
captures  the  steady  intensity  period  (not  always  an  easy  task 
for  mesoscale  models),  influenced  by  the  effect  of  ocean 
cooling  in  the  coupled  system.  The  ensemble  mean  and 
control  forecasts  are  similar  for  both  track  and  intensity, 
with  only  small  variations  for  different  forecast  times. 
Though  the  ensemble  mean  is  not  a  significant  improve¬ 
ment  over  the  control,  the  ensemble  does  provide  the 
possibility  for  uncertainty  information  about  the  track  and 
intensity. 

The  forecast  tracks  for  Hurricane  Ike  are  shown  in 
Fig.  7.  Both  coupled  and  uncoupled  ensemble  simulations 
are  performed  to  examine  the  impact  of  air-sea  coupling  on 
the  track  and  intensity  of  the  hurricane.  The  uncoupled 
simulation  is  conducted  in  an  identical  manner  to  the 
coupled  simulation  with  the  exception  that  there  is  no 
feedback  between  the  ocean  and  atmosphere.  Thus,  the 
analyzed  SST  computed  from  NCODA  is  held  fixed  in  time 
for  the  60-h  forecasts  for  the  14  uncoupled  ensemble 
simulations.  Generally,  the  forecast  tracks  for  the  time 
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period  0000  UTC  11  September  to  0000  UTC  13 
September  for  both  coupled  and  uncoupled  simulations 
are  similar,  with  the  majority  of  tracks  to  the  right  of  best 
track.  The  ensemble  statistics  for  the  coupled  simulation 
(Fig.  8)  clearly  indicate  that  this  coarse  resolution  simula¬ 
tion  is  unable  to  capture  the  intensity  of  the  observed 
hurricane  prior  to  landfall  (wind  speed  errors  ~2(M10  kts; 
SLP  errors  ~25-40  hPa  not  shown).  A  comparison  of  the 
control  and  the  ensemble  mean  indicates  that  the  ensemble 
mean  performs  no  worse  and  slightly  better  than  the  control  in 
predicting  track  and  intensity  (i.e.,  at  48  h,  mean  track  error  = 
105  vs.  125  nm  for  control).  The  use  of  the  coupled  model 
(versus  the  uncoupled)  has  a  general  positive  impact  on  track 
(Fig.  9).  The  coupled  ensemble  mean  track  error  is  generally 
smaller  than  the  uncoupled  track  error  (i.e.,  48-h  uncoupled 
track  error  =120  nm)  and  the  spread  is  generally  larger 
(48-h  coupled  ensemble  spread  —  53  vs.  45  nm  for 
uncoupled). 

Improvements  would  be  expected  at  higher  resolutions 
relative  to  this  coarse  resolution  (27  km).  It  should  be  noted 
that  high  resolution  (~5  km  or  less)  (Chen  et  al.  2007)  is 
typically  required  for  mesoscale  models  to  provide  any  skill 
in  forecasting  TC  intensity.  The  Hanna  and  Ike  cases 
provide  some  indication  that  the  ensemble  system  can  be  a 
useful  tool  to  improve  TC  track  and  intensity  forecasts, 
even  at  a  relatively  coarse  resolution. 

Figure  10  shows  the  SST  difference  between  0000  UTC 
11  September  and  0000  UTC  13  September  for  Hurricane 
Ike  for  the  TRMM  Microwave  Imager  satellite  observations 
and  the  coupled  ensemble  mean.  There  is  a  warm  core  eddy 


(WCE)  in  the  Gulf  of  Mexico  located  near  25°N,  87°  W  as 
indicated  by  the  sea  surface  height  contours  in  Fig.  10b. 
Observations  show  a  significant  SST  change  on  the  right 
side  of  the  storm  track  after  Ike  passes  the  WCE  with  SST 
decreased  by  up  to  ~4°C  over  48  h.  The  ensemble  mean 
shows  a  similar  location  for  the  SST  decrease  but  with  a 
smaller  magnitude  (~2°C).  The  area  of  cooling  in  the 
ensemble  mean  is  also  smaller.  The  difference  of  temper¬ 
ature  change  between  the  observations  and  the  ensemble 
mean  is  consistent  with  a  much  weaker  forecasted  storm 
intensity  as  compared  to  the  best  track  (Fig.  8a).  The  lack 
of  intensity  could  certainly  account  for  the  weaker  cooling, 
although  Hong  et  al.  (2000)  surmise  that  the  weaker  cooling 
could  be  related  to  modifications  in  the  pre-storm  ocean 
thermal  structure  as  a  result  of  changes  in  vertical  mixing. 
The  lack  of  cooling  indicates  that  the  upper  ocean 
stratification  is  not  strong  enough  to  represent  the  pre-storm 
conditions. 

During  the  life  cycle  of  a  TC,  the  region  along  the  TC  track 
should  also  be  a  region  of  large  ensemble  spread.  Figure  11 
shows  the  ensemble  spread  for  Hurricane  Ike  for  the 
24-h  forecast  valid  1200  UTC  11  September  of  surface  wind 
stress  and  surface  current.  The  ensemble  spread  of  the 
surface  wind  stress  is  maximum  in  the  relatively  large  region 
of  high  TC  wind  speeds  and  relatively  small  outside  the 
region.  Likewise,  the  spread  of  the  surface  current  is  a 
maximum  in  the  core  TC  region,  with  larger  values 
emanating  out  from  the  center  to  the  north-northeast  in  a 
clockwise  spiral.  The  region  of  large  spread  to  the  northwest 
of  the  Yucatan  Peninsula  is  believed  to  be  due  to  an  ocean 


Ensemble  spread 
24-h  forecast  valid  2008091 112 


Fig.  n  Ensemble  spread  for  the  29-member  Hurricane  Ike  case  for  the  24-h  forecast  valid  1200  UTC  11  September  of  a  surface  wind  stress 
(N  m2)  and  b  surface  current  (ms-1) 
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eddy  in  several  ensemble  members  and  is  believed  to  be 
spurious.  Spread-skill  relationships  of  atmosphere  and  ocean 
temperature  forecasts  at  various  levels  for  both  Hanna  and  Ike 
(not  shown)  further  support  the  under-dispersive  nature  of  the 
ensemble  system,  indicating  a  general  lack  of  spread  among 
ensemble  members.  This  is  a  common  feature  of  most 
ensemble  systems  (Houtekamer  and  Mitchell  2001),  and 
improved  techniques  are  being  developed  to  improve  the 
spread,  particularly  in  the  ocean,  and  will  be  presented  in  a 
follow-on  paper. 

6  Summary  and  conclusions 

The  development  and  testing  of  a  coupled  ocean-atmo¬ 
sphere  meso scale  ensemble  prediction  system  has  been 
described.  The  system  components  include  CO  AMPS  and 
NCOM  for  atmosphere  and  ocean  forecasting  and  NAVDAS 
and  NCODA  for  3D  variational  atmosphere  and  ocean  data 
assimilation.  The  components  are  efficiently  implemented 
under  the  Earth  System  Modeling  Framework  and  have  been 
rigorously  tested  for  a  variety  of  model  configurations  and 
case  studies.  The  Ensemble  Transform  technique  is  currently 
used  independently  for  the  ocean  (current,  temperature,  and 
salinity  variables)  and  atmosphere  (wind,  temperature,  and 
moisture  variables)  to  create  perturbations  for  the  coupled 
system.  Estimates  of  analysis  error  covariance  from  the  ocean 
and  atmosphere  3DVAR  assimilation  systems  are  used  to 
constrain  the  ET. 

Results  are  presented  for  tropical  cyclone  case  studies  of 
Hanna  and  Ike  (September  2008)  to  demonstrate  the 
baseline  capability  of  the  system.  A  suite  of  atmospheric 
perturbations  are  used  to  represent  the  uncertainty  due  to 
the  model  physical  parameterizations  (see  Table  1).  Gener¬ 
ally,  the  ensemble  mean  has  comparable  skill  as  compared 
to  the  deterministic  control  member  for  hurricane  track  and 
intensity  (even  at  a  relatively  coarse  27-km  resolution).  For 
Hanna,  the  ensemble  mean  forecast  track  errors  from  0  to 
48  h  are  well  within  the  average  errors  for  the  official 
forecasts  prior  to  landfall.  The  average  forecast  intensity 
error  among  members  is  ~10  kts  for  the  period  prior  to 
landfall,  while  the  official  NHC  forecasts  over-estimated 
Hanna  as  a  category-2  storm  upon  landfall.  The  average 
forecast  intensity  error  among  members  for  Ike  (-20-40  kts) 
is  greater  than  for  Hanna  as  both  the  coupled  and  uncoupled 
systems  did  not  initialize  the  storm  sufficiently  well  (initial 
errors  at  0000  UTC  1 1  September  of  -30  hPa).  Both  the  ocean 
and  atmosphere  components  of  the  ensemble  are  under- 
dispersive.  Efforts  are  underway  to  improve  the  spread, 
particularly  in  the  ocean. 

Future  plans  for  the  system  include  the  development  of 
ensemble  ocean-atmosphere  coupled  covariances  using 
innovative  localization  techniques.  This  approach  will 


enable  ensemble  4D-variational  (4DVAR)  coupled  data 
assimilation.  The  ensemble  4DVAR  system  will  provide 
improved  state  estimates  of  the  coupled  system  and 
increase  the  effectiveness  of  ensemble-based  targeted 
observation  approaches  that  are  currently  limited  by 
spurious  ensemble  correlations  and  by  neglected  error 
correlations  between  atmosphere  and  ocean  variables.  New 
case  studies  with  the  updated  system  are  proposed  for  both  TC 
and  non-TC  cases  to  investigate  and  quantify  the  skill  of  the 
coupled  ensemble  system  at  higher  spatial  resolutions.  It  is 
envisioned  that  the  new  updated  system  will  provide 
improved  mesoscale  ensemble  forecasts  for  use  in  probabilis¬ 
tic  products,  such  as  reliability  and  frequency  of  occurrence, 
and  in  risk  management  applications. 
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